Rational Design of Non-Covalent Imprinted Polymers Based on the Combination of Molecular Dynamics Simulation and Quantum Mechanics Calculations

Molecular imprinting is a promising approach for developing polymeric materials as artificial receptors. However, only a few types of molecularly imprinted polymers (MIPs) are commercially available, and most research on MIPS is still in the experimental phase. The significant limitation has been a challenge for screening imprinting systems, particularly for weak functional target molecules. Herein, a combined method of quantum mechanics (QM) computations and molecular dynamics (MD) simulations was employed to screen an appropriate 2,4-dichlorophenoxyacetic acid (2,4-D) imprinting system. QM calculations were performed using the Gaussian 09 software. MD simulations were conducted using the Gromacs2018.8 software suite. The QM computation results were consistent with those of the MD simulations. In the MD simulations, a realistic model of the ‘actual’ pre-polymerisation mixture was obtained by introducing numerous components in the simulations to thoroughly investigate all non-covalent interactions during imprinting. This study systematically examined MIP systems using computer simulations and established a theoretical prediction model for the affinity and selectivity of MIPs. The combined method of QM computations and MD simulations provides a robust foundation for the rational design of MIPs.


Introduction
Molecular imprinting technology (MIT) is a polymer-preparation-based synthesis process of artificial receptors for specific target molecules [1,2].Molecularly imprinted polymers (MIPs) prepared by MIT have the following significant advantages: low price, easy preparation, high selectivity and good reusability [3].Over the past few decades, MIT has successfully been employed in solid phase extraction [4][5][6], catalysis [7,8], drug delivery systems [9] and sensors [10,11].However, only a few types of MIPs are commercially available currently, and most research on MIP S is still in the experimental stage due to challenges in screening imprinting systems, particularly for weak functional target molecules.
Common screening methods for imprinting systems include experimental and ultraviolet (UV) spectral methods [12,13].Traditional experimental approaches primarily rely on numerous experiments and previous research.The accuracy of UV spectral methods is low [14,15].The introduction of computational chemistry and molecular simulation methods offers a new approach to screening imprinting systems [6,[15][16][17].Quantum mechanics (QM) computations can be employed to determine putative complexes among monomers and templates to obtain the optimal monomer with high selectivity and rebinding capacity to templates [18].Molecular dynamics (MD) simulations can simulate and predict the properties of an ensemble comprising non-covalent complexes formed between polymer components in a MIP pre-polymerisation mixture [14,19].There have been a few attempts for combining QM computations with MD simulations to predict the composition of pre-polymerisation mixtures for MIPs [20].
Herein, a combined method of QM computations and MD simulations was used to screen appropriate 2,4-dichlorophenoxyacetic acid (2,4-D) imprinting systems.This integrated screening method of imprinting systems leveraged the advantages of the two approaches to provide a more comprehensive description of the imprinting systems.QM computations were performed using the Gaussian 09 software, and MD simulations were performed using the Gromacs2018.8programme.Among various dynamics programmes, Gromacs provided the highest computational speed and performance and has improved consistently based on recent developments in Graphics Processing Unit (GPU) and Central Processing Unit (CPU) hardware [21,22].The steps involved in the MD simulations are illustrated in Figure 1.Based on the molecular structure of 2,4-D, 11 imprinting systems were examined, including two imprinting solvents (dimethyl sulfoxide (DMSO) and chloroform) and five functional monomers (acrylamide (AM), 2-(dimethylamino)ethyl methacrylate (DMAEMA), methacrylic acid (MAA), 2-trifluoromethacrylic acid (TFMAA) and 4-vinylpyridine (4-VP)).monomers and templates to obtain the optimal monomer with high selectivity and rebinding capacity to templates [18].Molecular dynamics (MD) simulations can simulate and predict the properties of an ensemble comprising non-covalent complexes formed between polymer components in a MIP pre-polymerisation mixture [14,19].There have been a few attempts for combining QM computations with MD simulations to predict the composition of pre-polymerisation mixtures for MIPs [20].Herein, a combined method of QM computations and MD simulations was used to screen appropriate 2,4-dichlorophenoxyacetic acid (2,4-D) imprinting systems.This integrated screening method of imprinting systems leveraged the advantages of the two approaches to provide a more comprehensive description of the imprinting systems.QM computations were performed using the Gaussian 09 software, and MD simulations were performed using the Gromacs2018.8programme.Among various dynamics programmes, Gromacs provided the highest computational speed and performance and has improved consistently based on recent developments in Graphics Processing Unit (GPU) and Central Processing Unit (CPU) hardware [21,22].The steps involved in the MD simulations are illustrated in Figure 1

QM Calculations for Characterising Template-Monomer Complexes
All geometric structures were computed using Gaussian 09 software [23] with the polarisable continuum model.The Becke 3-parameter-Lee-Yang-Parr (B3LYP) functional with a 6-311+G** basis set was used for optimisation based on Grimme's D3 dispersion correction [24,25].Frequency computations were conducted at the same level to determine whether a stationary point was a minimum at T = 298.15K. Interaction energy, ΔE, is computed using the following equation: where E(template-monomer complex) represents the energy of the complexes, E(template) denotes the energy of 2,4-D and E(monomer) represents the energy of functional monomers.

Computational Details 2.1. QM Calculations for Characterising Template-Monomer Complexes
All geometric structures were computed using Gaussian 09 software [23] with the polarisable continuum model.The Becke 3-parameter-Lee-Yang-Parr (B3LYP) functional with a 6-311+G** basis set was used for optimisation based on Grimme's D3 dispersion correction [24,25].Frequency computations were conducted at the same level to determine whether a stationary point was a minimum at T = 298.15K. Interaction energy, ∆E, is computed using the following equation: where E(template-monomer complex) represents the energy of the complexes, E(template) denotes the energy of 2,4-D and E(monomer) represents the energy of functional monomers.

MD Simulations for the Characterisation of the Types and Strengths of All Pre-Polymerization Interactions
All MD simulations were conducted using the Gromacs2018.8suite of programmes [21,22].The general amber force field [26] was used to establish and parameterise the molecular systems.The RESP20.5 charge method was used to assign the partial atomic charges.Multiwfn [27] used a highly efficient algorithm to evaluate the electrostatic potential employed in the analysis [28][29][30].The system composition is presented in Table 1.Table 1 shows the number of different molecules added during the MD simulations.Periodic boundary conditions were used in these simulations, and the simulation box's size was approximately 10 3 nm 3 .The order of simulation was as follows: 2,4-D (30 molecules) and a functional monomer (90-150 molecules) were added into a small box and solvated with 5597 molecules of chloroform/DMSO.This simulation condition approximately corresponded to 1 mmol of 2,4-D and 3-5 mmol of the functional monomer in 15 mL of solvent.Subsequently, the steepest descent approach was used to eliminate bad van der Waals contacts.Following these energy optimisation procedures, they were annealed for 2 ns under a constant volume at NVT, and temperature variations with simulation time were recorded (Figure S1).The boxes were then equilibrated at NPT (298.15K, 1 bar) for 3 ns.MD simulation production was performed for 5 ns at the NPT after equilibration of the systems.
For non-bonding interactions, a cut-off of 1.0 nm was used for all the prepared systems.The particle mesh Ewald simulation method was used to handle the electrostatic interactions in the systems.The cut-off simulation method was used to manage the van der Waals forces in these systems.During 5 ns of equilibration, a Berendsen barostat and velocity rescaling thermostat were used to control the pressure and temperature, respectively.During the 5 ns of production, the temperature and pressure were regulated using the Parrinello-Rahman barostat and velocity rescaling thermostat.
The trajectory and statistical analyses of the simulation data were performed using GROMACS' practical tools, and the results were visualised using visual molecular dynamics (VMD) [31].Hydrogen bond, radial distribution functions (RDFs) and spatial distribution functions (SDF) computations were primarily included in the analysis.RDFs were generated to measure the local concentrations of specific atom pairs at optimal interaction distances.The investigated atom pairs were described (Figure 2).Hydrogen bonds in Polymers 2024, 16, 2257 4 of 14 the trajectories were counted using a GROMACS command based on geometric standards.

O and N served as H-H bond acceptors and O-H and N-H served as H-H bond donors. The standard of H-H bonds is introduced below. The angle of the H-H donor-acceptor was less than 30
• and the distance between the donor and the acceptor was less than 3.5 Å (Figure S2).The SDF was generated using the gmx spatial tool in the GROMACS and VMD packages.First, the trajectory was processed using the gmx spatial tool in GROMACS to create a cube file.To obtain meaningful SDFs, the solute molecules must be centred within the box for the entire trajectory, and the molecules for which SDF computations are required should be grouped and defined.Subsequently, the SDF is visualised using the VMD package.
The trajectory and statistical analyses of the simulation data were performed using GROMACS' practical tools, and the results were visualised using visual molecular dynamics (VMD) [31].Hydrogen bond, radial distribution functions (RDFs) and spatial distribution functions (SDF) computations were primarily included in the analysis.RDFs were generated to measure the local concentrations of specific atom pairs at optimal interaction distances.The investigated atom pairs were described (Figure 2).Hydrogen bonds in the trajectories were counted using a GROMACS command based on geometric standards.O and N served as H-H bond acceptors and O-H and N-H served as H-H bond donors.The standard of H-H bonds is introduced below.The angle of the H-H donor-acceptor was less than 30° and the distance between the donor and the acceptor was less than 3.5 Å (Figure S2).The SDF was generated using the gmx spatial tool in the GROMACS and VMD packages.First, the trajectory was processed using the gmx spatial tool in GROMACS to create a cube file.To obtain meaningful SDFs, the solute molecules must be centred within the box for the entire trajectory, and the molecules for which SDF computations are required should be grouped and defined.Subsequently, the SDF is visualised using the VMD package.

QM Calculations for Characterising Template-Monomer Complexes
The optimised conformations of 2,4-D and the five functional monomers were obtained (Figure S3).The optimised conformations of the complexes are illustrated in Figure 3.The interaction energies of the 2,4-D-functional monomer complexes in different solvents are presented in Tables 2 and 3.As shown in Table 2, the stability of the complexes in chloroform decreased in the following order: MAA (−63.68 kJ/mol) > TFMAA (−63.00 kJ/mol) > AM (−61.82 kJ/mol) > 4-VP (−58.54 kJ/mol) > DMAEMA (−39.31kJ/mol).As shown in Table 3, the stability order of the complexes in DMSO was as follows: TFMAA (−57.84 kJ/mol) > AM (−56.17 kJ/mol) > 4-VP (−55.56 kJ/mol).When the porogen was changed from CHCl3 to the more polar DMSO, the binding stability between 2,4-D and the functional monomers decreased.
To understand the recognition method between functional monomers and 2,4-D in chloroform, the molecular electrostatic potential (MESP) mapped on molecular van der Waal (vdW) surfaces of 2,4-D and the active monomers was calculated (Figure 4).The variation in the mapped MESP is illustrated using the following coding scheme: red represents the most positive potential regions, blue indicates the most negative potential regions and white area signifies potential halfway between the two extremes of blue and red.Van der Waals surface penetration was visible in the region where hydrogen bonds were formed, and red and blue regions penetrated each other (Figure 4).This indicated the complementary characteristics of electrostatic potential and the essence of the electrostatic attraction interaction of the H-H bond.
Polymers 2024, 16, 2257 6 of 14 resents the most positive potential regions, blue indicates the most negative potential gions and white area signifies potential halfway between the two extremes of blue a red.Van der Waals surface penetration was visible in the region where hydrogen bon were formed, and red and blue regions penetrated each other (Figure 4).This indica the complementary characteristics of electrostatic potential and the essence of the elect static attraction interaction of the H-H bond.

Validity of the Force Field and Thermostat
To confirm the validity of the force field, density was considered as a key prope The comparison of experimental and simulated densities at 298.15 K is shown in Tabl At this temperature, the computed densities were consistent with the experimental val and the error was ~5%, revealing that the constructed molecular force field could rep duce the experimental density.The force field was verified through property compu tions.In addition, temperature fluctuations within 5 ns of production were recorded (F ure S4).Table 4. Calculated density of organic small molecules under force fields (unit: 1 g/cm 3 ).

MD Simulations for Characterising the Types and Strengths of All Pre-Polymerisation Interactions 3.2.1. Validity of the Force Field and Thermostat
To confirm the validity of the force field, density was considered as a key property.The comparison of experimental and simulated densities at 298.15 K is shown in Table 4.At this temperature, the computed densities were consistent with the experimental values and the error was ~5%, revealing that the constructed molecular force field could reproduce the experimental density.The force field was verified through property computations.In addition, temperature fluctuations within 5 ns of production were recorded (Figure S4).DMSO, a universal solvent, is often selected as an imprinting solvent.The configurations of the last frames in all imprinted systems (P1-P3) with DMSO as a solvent are shown in Figure 5. Numerous DMSO molecules were wrapped with 2,4-D and active monomers.The DMSO molecules in Figure 5 were temporarily masked to enable more obvious observations of the interaction between 2,4-D and the active monomer (Figure S5).The functional monomer and 2,4-D's hydrogen bond interactions were almost completely disrupted by the addition of DMSO, a protic polar porogen (Figure S5).tions of the last frames in all imprinted systems (P1-P3) with DMSO as a solvent are shown in Figure 5. Numerous DMSO molecules were wrapped with 2,4-D and active monomers.The DMSO molecules in Figure 5 were temporarily masked to enable more obvious observations of the interaction between 2,4-D and the active monomer (Figure S5).The functional monomer and 2,4-D's hydrogen bond interactions were almost completely disrupted by the addition of DMSO, a protic polar porogen (Figure S5).To evaluate the interactions between the functional monomer and 2,4-D for each imprinting system using DMSO as a solvent, the total non-bonding interaction was divided into van der Waals (EvdW) and electrostatic (Eelec) (Table 5).In each imprinting system with DMSO, the binding energy between 2,4-D and various functional monomers was small.
RDFs were used to determine the average radial distribution of functional monomers around 2,4-D based on the trajectories.The RDFs indicated that the various functional monomers had a low distribution density in the vicinity (<3.0 Å) of 2,4-D and were almost in the average distribution state in the system (Figure 6).Furthermore, no hydrogen bonds were formed between the functional monomers and 2,4-D (Figure S6).The 2,4-D was not well combined with functional monomers, displaying the free dispersion state, because the polar solvent molecule, DMSO, significantly interfered with the electrostatic force, such as hydrogen bonding and electrostatic interaction.
RDFs were used to determine the average radial distribution of functional mon around 2,4-D based on the trajectories.The RDFs indicated that the various func monomers had a low distribution density in the vicinity (<3.0 Å) of 2,4-D and were a in the average distribution state in the system (Figure 6).Furthermore, no hydrogen were formed between the functional monomers and 2,4-D (Figure S6).The 2,4-D w well combined with functional monomers, displaying the free dispersion state, be the polar solvent molecule, DMSO, significantly interfered with the electrostatic such as hydrogen bonding and electrostatic interaction.

Imprinting Systems Using Chloroform as the Solvent
Chloroform, a non-polar solvent, is widely employed in the production of non-covalent MIPs.The conformations of the last frames of all imprinted systems (P4-P8) using chloroform as a solvent are shown in Figure 7.To better observe the interactions between 2,4-D and the functional monomers, the chloroform molecules shown in Figure 7 were temporarily masked (Figure S7).The combination of functional monomers and 2,4-D was significantly affected after the solvent was changed from DMSO to CHCl 3 (Figure S7).Throughout the simulation of 5 ns, AM, MAA and TFMAA exhibited stronger and identical interactions with 2,4-D.DMAEMA was involved in lower levels of binding with 2,4-D and the 4-VP-2,4-D interaction did not exist.The differences in behaviour can be attributed to the relative strengths of the functional groups of the five functional monomers.The primary modes of interaction between 2,4-D and the different functional monomers were further assessed (Figure S8).The binding between the functional monomers and 2,4-D was hydrogen bonding.The 2,4-D, AM, MAA and TFMAA demonstrated a higher degree of dimerisation in the non-polar solvent CHCl 3 because amide nitrogen and carboxyl groups could form two hydrogen bonds.
The interaction energies between the various functional monomers and 2,4-D in the imprinting systems using chloroform as the solvent are presented in Table 6.To compute the interaction energies, the total non-bonding interaction energies between the functional monomers and 2,4-D were divided into EvdW and Eelec.The non-bonding interaction (−1221.70kJ/mol) between AM and 2,4-D was significantly stronger than that between 2,4-D and the other functional monomers.
Polymers 2024, 16, 2257 9 of 14 tributed to the relative strengths of the functional groups of the five functional monomers.The primary modes of interaction between 2,4-D and the different functional monomers were further assessed (Figure S8).The binding between the functional monomers and 2,4-D was hydrogen bonding.The 2,4-D, AM, MAA and TFMAA demonstrated a higher degree of dimerisation in the non-polar solvent CHCl3 because amide nitrogen and carboxyl groups could form two hydrogen bonds.The interaction energies between the various functional monomers and 2,4-D in the imprinting systems using chloroform as the solvent are presented in Table 6.To compute the interaction energies, the total non-bonding interaction energies between the functional monomers and 2,4-D were divided into EvdW and Eelec.The non-bonding interaction (−1221.70kJ/mol) between AM and 2,4-D was significantly stronger than that between 2,4-D and the other functional monomers.RDFs were used to determine the mean radial distribution of the functional monomers around 2,4-D from the resulting trajectories.RDFs indicated that the probability of AM being close (<3.0Å) to 2,4-D was higher than the mean probability (Figure 8).The number of H-H bonds formed between the functional monomers and 2,4-D is shown in Figure 9. AM formed a significantly higher number of H-H bonds with 2,4-D (~45) than those formed by the other functional monomers.Meanwhile, very few hydrogen bonds (<5) were formed between 2,4-D and the basic functional monomers.mers around 2,4-D from the resulting trajectories.RDFs indicated that the probabil AM being close (<3.0Å) to 2,4-D was higher than the mean probability (Figure 8 number of H-H bonds formed between the functional monomers and 2,4-D is sho Figure 9. AM formed a significantly higher number of H-H bonds with 2,4-D (~45 those formed by the other functional monomers.Meanwhile, very few hydrogen b (<5) were formed between 2,4-D and the basic functional monomers.To examine the spatial distribution of AM, the SDF of AM around 2,4-D was statistically analysed (Figure S9).The density grid representations revealed that the carboxyl group was the primary binding site of 2,4-D in hydrogen bond formation.

Effect of the Ratio of the Template to Functional Monomer
The molar ratio of the template to functional monomer significantly influenced the formation of the template-functional monomer complexes.This study further investigated the effects of the ratio of 2,4-D to AM based on the P8 imprinting system.Three ratios (namely 1:3, 1:4 and 1:5) were selected.The conformations of the last frames of the imprinting systems (P8-P10) are shown in Figure 10.For an enhanced observation of the To examine the spatial distribution of AM, the SDF of AM around 2,4-D was statistically analysed (Figure S9).The density grid representations revealed that the carboxyl group was the primary binding site of 2,4-D in hydrogen bond formation.The molar ratio of the template to functional monomer significantly influenced the formation of the template-functional monomer complexes.This study further investigated the effects of the ratio of 2,4-D to AM based on the P8 imprinting system.Three ratios (namely 1:3, 1:4 and 1:5) were selected.The conformations of the last frames of the imprinting systems (P8-P10) are shown in Figure 10.For an enhanced observation of the interactions between 2,4-D and AM, the chloroform molecules in Figure 10 were temporarily masked (Figure S10).As the AM content increased, the amount of free AM in the systems also increased, leading to the formation of non-specific adsorption sites.The number of hydrogen bonds formed between 2,4-D and AM in the imprinting systems (P8-P10) is shown in Figure 11.With the increase in AM content, the number of hydrogen bonds formed between AM and 2,4-D did not increase significantly, but increased among AM and AM (Figure 11).molecules (black); TFMAA molecules (green); MAA molecules (green); AM molecules (cyan)).To better observe the binding between 2,4-D and the functional monomers, chloroform molecules were temporarily masked.Figure S8.Main mode of the complexes in the last frame of the imprinting systems P4-P8.Figure S9.Three-dimensional grid-density representations of the imprinting system P8.Densities describe the probability and associated lifetimes for finding AM around 2,4-D.Higher contour values for AM are presented in blue.Figure S10.Conformations of the last frames of the imprinting systems P8-P10.(2,4-D molecules (red); AM molecules (cyan)).To better observe the binding between 2,4-D and AM, chloroform molecules were temporarily masked.

Figure 1 .
Figure 1.Schematic of the steps involved in a typical molecular dynamics simulation for a prepolymerisation mixture illustrated for a system based on a functional monomer (AM) and 2,4-D in chloroform: (A) Components added to the design.(B) A 5 ns production phase and the pre-polymerisation mixture after equilibration.(C) Statistical analysis results of the different pre-polymerisation mixtures produced.

Figure 1 .
Figure 1.Schematic of the steps involved in a typical molecular dynamics simulation for a prepolymerisation mixture illustrated for a system based on a functional monomer (AM) and 2,4-D in chloroform: (A) Components added to the design.(B) A 5 ns production phase and the prepolymerisation mixture after equilibration.(C) Statistical analysis results of the different prepolymerisation mixtures produced.

Figure 6 .
Figure 6.RDFs displaying probabilities of obtaining the atomic densities of the functional mers at various separation distances from the 2,4-D functional groups in the imprinting syste P3.

Figure 6 .
Figure 6.RDFs displaying probabilities of obtaining the atomic densities of the functional monomers at various separation distances from the 2,4-D functional groups in the imprinting systems P1-P3.

Figure 8 .
Figure 8. RDFs displaying probabilities of obtaining the atomic densities of the functional mers for various separation distances from the 2,4-D functional groups in the imprinting sy P4-P8.

Figure 8 .
Figure 8. RDFs displaying probabilities of obtaining the atomic densities of the functional monomers for various separation distances from the 2,4-D functional groups in the imprinting systems P4-P8.Polymers 2024, 16, x FOR PEER REVIEW 11 of 15

Figure 9 .
Figure 9. Number of H-H bonds produced between the functional monomers and 2,4-D in the imprinting systems P4-P8 with chloroform as the solvent.

Figure 9 .
Figure 9. Number of H-H bonds produced between the functional monomers and 2,4-D in the imprinting systems P4-P8 with chloroform as the solvent.

4 .
Effect of the Ratio of the Template to Functional Monomer

Table 1 .
System design and molecular dynamics methodology used.

Table 4 .
Calculated density of organic small molecules under force fields (unit: 1 g/cm 3 ).

Table 5 .
Interaction energies (kJ/mol) of 2,4-D and different functional monomers in DMSO based on MD simulations.Calculated at 298.15 K and a pressure of 1 bar.

Table 6 .
Interaction energies (kJ/mol) of 2,4-D and the different functional monomers in chloroform based on MD simulations.Calculated at 298.15 K and a pressure of 1 bar.

Table 6 .
Interaction energies (kJ/mol) of 2,4-D and the different functional monomers in chloroform based on MD simulations.Calculated at 298.15 K and a pressure of 1 bar.